#include "operator.h"
void X_Lambda::init(int hw) {
  this->hbar_omega = hw;
  cout << "X_Lambda hw = " << this->hbar_omega << endl;
  int N_max = 16;
  frac_vec.resize(N_max + 1);
  frac_vec[0] = 1;
  for (int i = 1; i <= N_max; i++) {
    frac_vec[i] = frac_vec[i - 1] * i;
    // cout<<" i = "<<i<<" \t frac : "<<frac_vec[i]<<endl;
  }
}

double X_Lambda::cal(int n_p, int l_p, int n, int l, int Lambda) {
  double res = 0.0;
  int tau_2 = (l_p - l + Lambda);
  int tau_p_2 = (l - l_p + Lambda);
  if (tau_p_2 < 0 || tau_2 < 0 || tau_p_2 % 2 != 0 || tau_2 % 2 != 0) {
    return res;
    cout << "Wrong happened @  X_Lambda::cal " << endl;
    exit(0);
  }
  int tau_p = tau_p_2 / 2;
  int tau = tau_2 / 2;
  double t = 0.5 * (l_p + l + Lambda + 1);
  //
  int sigm_min_t = std::max(n - tau, n_p - tau_p);
  int sigm_min = std::max(sigm_min_t, 0);
  int sigm_max = std::min(n, n_p);
  if (sigm_max < sigm_min) {
    return res;
  }
  //
  double C1_t1 = gsl_sf_gamma(n + 1.0) * gsl_sf_gamma(n_p + 1.0);
  double n_d = n + 1 + t - tau;
  double n_p_d = n_p + 1 + t - tau_p;
  // if (n_b_d < 0 || n_a_d < 0) {
  //   cout << "Wrong happened @  X_Lambda::cal 2" << endl;
  //   exit(0);
  // }
  double C1_t2 = gsl_sf_gamma(n_d) * gsl_sf_gamma(n_p_d);
  double C1 = std::sqrt(C1_t1 / C1_t2);
  //
  double C2 = frac_vec[tau_p] * frac_vec[tau];
  //
  double C3 = 0.0;
  for (int s = sigm_min; s <= sigm_max; s++) {
    double C3_1 = gsl_sf_gamma(t + s + 1.0);
    double C3_2 = frac_vec[s] * frac_vec[n - s] * frac_vec[n_p - s] *
                  frac_vec[s + tau - n] * frac_vec[s + tau_p - n_p];
    C3 += C3_1 / C3_2;
  }
  res = C1 * C2 * C3 * hbar_omega;
  return res;
}
